arXiv:astro-ph/0211026vl 1 Nov 2002 


Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 2 February 2008 (MN lAT^^X style file v2.2) 


The impact of mass loss on star cluster formation. II. 
Numerical N-body integration & further applications 

C. M. Boily^’^ & P. Kroupa^ 

^Astronomisches Rechen-Institut, Monchhofstrasse 12-14 Heidelberg, D-69120 Germany 

^Institut fiir Theoretische Physik und Astrophysik der Universitdt Kiel, D-24098 Kiel, Germany 

^Present address: Observatoire astronomique de Strasbourg, 11 rue de I’universite, 67000 Strasbourg, France 


2 February 2008 


ABSTRACT 

We subject to an N-body numerical investigation our analysis of Paper I on the sur¬ 
vival of stellar clusters undergoing rapid mass loss. We compare analytical tracks of 
bound mass-fraction vs star formation efficiency e to those obtained with N-body in¬ 
tegration. We use these to argue that stellar clusters must develop massive cores of 
high-binding energy if they are to remain bound despite a star formation efficiency as 
low as 30% or lower suggested by observations. The average local virial ratio (ct^/|^|) 
is introduced to classify bound clusters as function of their critical e for dissolution. 
Clusters dissolving at lower e achieve the lowest ratio. We applied this classification 
parameter successfully to Michie-King and Hernquist-type distribution functions. The 
Plummer sphere is exceptional in that it defies this and other classification parameters 
we tried. The reasons for the discrepancy include less effective energy redistribution 
during the expansion phase for this case. 
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1 INTRODUCTION 

This article is the second of two addressing the question of the survival of bound star clusters which suffer mass loss through 
gas expulsion. Rapid gas removal may unbound a star cluster if the star formation efficiency e falls below a critical value, 
where 


_ M. M. 

^ M M* -b Mgas ^ ’ 

is the ratio of stellar mass M* to total system mass M (cf. Hills 1980; Adams 2000; Clarke et al. 2000). The SFE t < 0.4 at 
birth derived observationally on the scales of young clusters (Lada 1999). Using the virial theorem. Hills (1980) found that 
the gas-free cluster would settle to a new configuration of radius i?* known in terms of the initial radius R from 

R. _ I _ 1 e 

R^2M*-iM^2e-i ^ ’ 

and hence > oo if e < by and large the value of reference in this field (e.g., Clarke et al. 2000). In Boily & Kroupa 
(2002, hereafter Paper I) we have shown that the value e for which R* ^ oo is function of the equilibrium stellar distribution 
function (DF) of the cluster at formation. The fraction of stars lost is obtained from the phase-space DF F{r, v) through the 
integrals (Paper I) 


1-Ae = 


rve rve 
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/ F{r,v) d^r du / f{v) 
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where f{v)A^ v = F{r, ?;)d®r d^u is the stellar velocity DF at birth, and Ae the bound fraction of stars following gas expulsion. 
Iterating on (^ to account for escaping stars provides a self-consistent selection criterion based on the local escape velocity, 
Ue,* (see Paper I for further details). 
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The procedure outlined here does not account for the dynamical evolution triggered by the loss of gas, when, on the 
whole, the residual gas mass {Mgas above) represents a non-negligible fraction of the total system mass M. In reality stars 
must orbit in space before leaving the system, and they may exchange gravitational energy with the background potential in 
the process (cf Binney & Tremaine 1987 = BT+87 for a discussion), when the cluster expands as a result of the shallower 
potential well felt by the stars. The purpose of this article is to establish the validity of (^) once dynamical evolution is fully 
taken into account. In this respect we continue in the lineage of Lada, Margulis & Deardorn (1984) but with the benefice of 
much larger-N simulations than was possible then (tens of thousand compared with hundreds). This allows us to work in the 
strictly collision-less regime, with much reduced noise in the potential. Because we can bridge over to Paper I, we will focus 
on a few calculations while exploring a wide range of DF. 

We approach the dynamics through N-body numerical integration as described in the following section. We compare 
results for Plummer spheres with the analytic result of Paper I. We then apply (^ to truncated isothermal (King) DF, as well 
as centrally peaked Hernquist-type models in subsequent sections. We derive limits from the analytic treatment applicable to 
star clusters, namely that embedded clusters with an SFE as low as 30% require a significant fraction of stellar orbits on low 
velocity (high binding energy) in order to withstand rapid gas evacuation, concurring with the detection of high-mass stars 
preferentially in dense cluster cores (Testi et al. 1997, 1999). We explore alternative mechanisms to form bound open clusters 
in the closing section of the paper. 


2 DETAILS OF THE NUMERICAL SETUP 

We performed N-body calculations with the multi-grid FFT code SUPERBOX (Fellhauer et al. 2000) and 50,000-particle models. 
We have chosen the scale of N-body models such that the total mass M = 1, G = 2, and cutoff radius R = 1.25 in all cases. The 
grid resolution and the selection of time-step was done such that the equilibria remained stable, in the sense that concentric 
10% mass radii remain constant upon integrating for 5 model units of time (from here on we write m.u. for model units), 
corresponding to ~ 6.73 ter where the crossing time ter is defined from the mean density p = , 


t 


cr 



~ 0.75 m.u. . 


( 4 ) 


Typical grid resolutions of from 0.005 to 0.01 in the central region proved adequate for the range of models considered. The 
same remark applies to all simulations discussed here. 

The nature of the problem at hand means we do not have to treat the gas dynamics as such: instead we note that the 
virial ratio 


2T M{v'^) _ 

^ ^ ~ GM'^/R MIR ’ 


( 5 ) 


where T and W are the kinetic and gravitational energies, respectively (Q = 1 defines virial equilibrium), relates velocity 
dispersion, mass and radius uniquely. Since gas dynamics is only taken into account insofar as it sets the stellar velocity 
dispersion, we use (^ to multiply equilibrium particle velocities by the factor M/M* = + Mgas/M* = as defined 

in (^). As a result Q > 1 when we replace M Mi, in (^: the super-virial system is dynamically identical to the problem 
posed. 

The excess kinetic energy leads to rapid expansion of the star cluster. We followed the evolution of the mass distribution 
in two ways. First we monitored concentric spherical shells, each enclosing 10% of the total stellar mass. Second, we counted 
the number of stars that remained within a volume of linear size twice as large as the initial system radius, or 2 x 1.25 = 2.5 
in model units. We found both approaches yield similar results. 

In the first instance the fraction of bound stars at the end of the calculation is found to 5% accuracy by plotting the 
Lagrange radii as function of time and determining which ones have turned around and settled to finite values. In the second 
instance, a star whose radial position shifted from R to 2R would have suffered a drop in potential from <()* = GM*/i7 to 
(pi, = GMi,/2R if no stars were lost, which corresponds to the mean binding energy per unit mass of the initial cluster of 
radius R. Thus counting as unbound stars orbiting outside 2R corresponds to removing stars whose binding energy is not 
larger than the mean value derived from the gravity of the stars alone. Both methods over-estimate the fraction of stars lost 
since they impose their own selection of bound orbits by binding energy: consequently, the numerical results provide lower 
limits to the fraction of bound stars at given SFE. We allowed stars more time for escape by doubling the integration time of 
a few calculations. In all cases our results did not vary by significant amounts, so we adopted as standard the setup described 
above. 
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3 PLUMMER MODELS 
3.1 Analytic solution 

We recall the analytic result of Paper I for Plummer spheres before discussing the time-evolution and bound fraction obtained 
from N-body runs. The isotropic DF of the n = 5 polytropic Plummer sphere is F{E) oc (—= (—Inserting 
this in (|^ leads to the parametric solution 

Ae(x) = — ^sin“^ — p{x) , (6) 

where p{x) = 1 — 1210/105 x + 2104/105 — 1488/105 x^ + 384/105 x‘^ and the parameter x 

x = \,e- a; C [0,1] . (7) 

The full solution Ae(a;),e(a:) follows from (^, known function of x, i.e. e{x) = x/\^{x) . The solution curve e{x) reaches a 
minimum at a; « 0.2252 when e ~ 0.442. Solutions with smaller e are impossible, which we interpret as indicating total cluster 
dissolution. Note that the solution A(a:) is independent of radius, which means that the same fraction of stars are escaping at 
each radius. 


3.2 N-body calculations 

The N-body realisation of a Plummer model is complete once we specify a truncation radius and a Plummer length Rp. We 
set Rp = 0.10 so that the numerical models contain 99% of the total mass of the model within r — 1.25. Note that the mean 
density of a Plummer sphere is (Rp/r)^ ~ 10~®x the central density so that ter corresponds to about thirty central crossing 
times: consequently the timesteps of integration was set to < 1/30 m.u. to resolve the dynamics at the centre. 

The results of runs with e ranging from 0.4 to 1 are listed in Table Ih (see also Fig. 2 of Paper I). We find the general 
trend and quantities in good agreement with the results derived from (™. Specihcally, the fraction of bound stars drops 
very rapidly around e = 0.44; furthermore, we find no indication that stars remain bound for the case where e = 0.40. The 
N-body calculation for e = 0.45 gives a bound fraction of ~ 44%. The critical SFE from our numerical models is therefore 
42.5 ± 2.5 %. This is in agreement with the findings by others for similar initial conditions (Lada, Margulis & Deardorn 1984; 
Goodwin 1997; Geyer & Burkert 2001). We may draw some confidence from the agreement between N-body models and the use 
of different algorithms for integration: for example, Geyer & Burkert used a TREE algorithm for their SPH-Nbody calculations. 

Table ^ lists the ratios between radii of hve concentric spherical shells taken when the system has settled to a new 
equilibrium, and their respective initial radius. Each shell therefore encloses the same Lagrangian mass as it did initially, i.e. 
the same fraction of the initial total stellar mass; we therefore refer loosely to the shells as ‘Lagrangian radii’. The average 
of these five ratios is given in the column noted E*/Eo and compared with equation (^. Values derived from (^ typically 
over-estimate the mean cluster radius, however note that the disagreement with the numerical N-body calculations becomes 
severe only for low SFE. The ratios of Lagrange radii help assess the morphological evolution of the cluster in the following 
way: If the initial Plummer profile after gas expulsion and dynamical relaxation were to settle to an expanded (homologous) 
version of itself, the ratios between the Lagrange radii at the end of the calculations and their initial values would give the 
same expansion factor for every Lagrange radius. At the end of the N-body calculations, in equilibrium, we Hnd this to hold 
true in the inner region of the system and for an SEE > 70%. We verified that these fluctuations do not originate from 
numerical resolution problems: test calculations of equilibrium models indicated displacements of mass shells of a few percent 
at most. At fixed SFE, the variations between individual ratios do not exceed 10% up to the 70% mass shell. As the SFE is 
reduced, the trend of near-constant expansion factor is gradually lost, such that when the SFE « 60% or lower the expansion 
factor increases monotonically with Lagrangian radii. The outermost mass shell expands the most in all the cases, until the 
SFE becomes so low that much of the initial mass leaves the simulation volume. Therefore, the new equilibria established by 
the star clusters differ signihcantly from the initial Plummer profile, in the sense that its core expands less than the outer 
envelope, and is not, as a result, an homologous map of the initial mass profile. (Figure ^[b,c] makes this point graphically, 
as we discuss below.) 

Equation (^) implicitly assumes that the number of stars with energy E > 0 is independent of details of the time-evolution. 
Comparisons of bound fractions with N-body calculations partly justifies this. To see why the iterated selection criterion (^ 
works, we counted the number of stars with positive binding energy E as function of time, as well as the mean radius of two 
subsets of stars: those with negative energy (E < 0), and all the stars within the simulation volume. 

The results are displayed on Fig. ^(a) for a Plummer model with an SFE = 50% (or, Q = 2). A key feature on this hgure is 
the quick leveling off of the fraction of E > 0 stars at t « 0.2 m.u. (shown as solid line on the figure), while both radii expand 
until eventually they reach a maximum at f « 0.75 (dashed lines). The number of E > 0 stars increases by a factor « 2.5 in 
the course of evolution, indicating a fair amount of energy exchanges soon takes place. A similar factor can be computed from 
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(^) and (^). The instantaneous fraction of unbound stars when the gas mass fraction is removed is 1 — Ae(e = 1/2) = 0.0877 
from (P), while the net fraction of escapers after iteration is 1 — Xe{x) = 0.15 with e(x = 0.425) = 1/2 from (^). The ratio 
0.15/0.0877 « 1.7 < 2.5 obtained from the N-body calculation. Nevertheless, the analytic iterative scheme encapsulates the 
basic dynamics at work during the expansion phase. 

The time-evolution of the system is characterised by the two mean radii (Fig. ^[a], dashed lines). At early times, both 
radii peak simultaneously, however they reach different values as might be expected from the selection criterion. This suggests 
that stars on the whole follow a global radial expansion, and hence few E < 0 stars have turned around before the expansion 
phase is over. 

The distribution of binding energy provides insight into the speed at which Plummer spheres evolve. Fig. |^(b,c) shows 
two snapshots of the energy distribution at times t = 1/2 m.u. and 2, respectively. We have graphed the particles’ energy 
versus their initial energy (after gas expulsion) to highlight evolution. The dash triangle isolates a number of stars which had 
negative binding energy at t = 0 but now have positive energy due to evolution. Note that none of the stars that had positive 
energy initially acquires negative binding energy at later times: indeed their energy at any times parallels the diagonal set 
by the initial conditions, as would be the case for particles cruising away from the cluster centre and not participating in the 
evolution of the central region. On Fig. ic) a clump forms with high negative binding energy, demonstrating that the cluster 
is still in its early stages of relaxation at that time. Yet the number of stars with E > 0 remains practically unchanged at 
around 11500 from then onwards (solid line. Fig. |^[a]). The non-homologous character of the dynamics is confirmed when we 
consider the changes in shape of the distribution of particles at different times. 

The numerical N-body experiments with Plummer spheres show that cluster dissolution takes place at SFE « 42.5±2.5%, 
in good quantitative agreement with (^. In Paper I, we showed that polytropes of index n > 5 would dissolve for lower SFE. 
However such polytropes are infinite in mass and radius, and hence they are of limited use. One exception is the truncated 
isothermal sphere (Michie-King models, of index n ^ oo) which have become benchmarks in globular cluster studies. 


4 MICHIE-KING MODELS & AN ALTERNATIVE SCHEME 
4.1 Analytic solution 

The luminosity profiles of globular star clusters are well fitted by Michie-King models (Michie 1963, King 1966; see Meylan 
& Heggie 1997 for a recent review). The truncated isothermal DF of Michie-King models takes the form (BT±87, §4.4) 


F{E) 


pi(27ra2)-®/2 



\<f> - 

n-2 




( 8 ) 


for V < ^J—2{(j> — 4>Vt]) and E = 0 otherwise. All Michie-King models are truncated at a finite radius, rt, where the density 
drops to zero. They are differentiated from one another by the central ratio 


which is a free parameter. As a function of radius, the mass density p flattens out near the centre in all cases. The volume of 
near-constant density bounded by the core radius Vo 


Vo = (9cr^/47rGp[4'/cr^]) (9) 

shrinks as tk/o-^ increases. We find on inserting (|^ in the integral (|^ 



where />* is the gravitational potential of the stars alone and we have used dimensionless variables 


(f){r) = </(r)/4'; v{r) = v(r)/a . 

Note that 1 — Xe{<j>) depends on the local potential, and hence it is a function of radius. To compute the net fraction of bound 
stars, we must therefore re-compute the potential numerically for each evaluation of Ae in (^). This poses no problem since 
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4>{r oo) ^ 0, and the density is known at each step (though it is no longer a King-Michie profile). The integral ( [lo[ ) was 
evaluated from standard algorithms (Press et al. 1992, p. 213) while the potential was computed by mid-point averages on a 
mesh of from 3000 to 6000 points, which proved sufficiently accurate for convergence of the results and in particular to identify 
the critical SFE for dissolution. Convergence was defined as variations 5(1 — Ae) averaged over the mesh smaller than one part 
in a million for two successive iterations. Results are shown on Fig. ^ for an T/cr^ = 9 King model, where we graphed the 
density (top panel) and bound fraction Ae as function of radius for three different values of the SFE, e. The radial dependence 
of Ae is accentuated for lower values of the SFE, however it always decreases monotonically with radius until it levels off again 
near r = rt = 1.25. For an SFE = 45%, we find Ae(0) ~ 0.90 at the centre, and Xe{rt) ~ 0.40 near the edge. The profile Ae(?') 
remains flat inside r = ~ 0.013 for this model, independently of the SFE. The core region of King models is therefore more 

robust against dissolution. By contrast, the solution for the Plummer model remains flat at all radii (Fig. la]). 


4.2 N-body calculations 

We compared the analytic prediction ( 0 ) with numerical N-body calculations for an = 9 model. Our choice of param¬ 

eter is motivated by the fact that the binding energy of Michie-King models peaks at approximately this value (all models of 
constant mass scaled to the same truncation radius rt), which should, therefore, provide the best hope to preserve a bound 
cluster despite low SFE. 

Numerical N-body calculations were done similarly to the case of Plummer spheres, taking care to adjust the grid 
resolution to achieve the same mass resolution per mesh. This way we did not bias the calculations against the more centrally 
peaked Michie-King models. The scales of time, mass and length are the same in both the cases. 

Fig. |^(a) graphs the bound mass fraction as function of the SFE. The dash is the analytic result ( 0 ), while the open 
circles show the N-body results. For an SFE = 50%, we integrated one more N-body model up to t = 25 m.u. to allow for a 
possible larger number of stars to escape, which would bring down the bound fraction. The two open circles seen at SFE = 
50% on Fig. ^(a) differ by a small amount however, at Ae = M*/M* = 0.64 (t = 5) and 0.59 [t = 25), respectively; this gives 
indirect confirmation of reduced evolution for t > 5 and boosts confidence in the analysis. The dashed curve does not give 
a good fit to the N-body data for the full range of SFE. We may understand the reason for this discrepancy Iw comparing 
the evolution of stars with negative energy to the cluster as a whole, as done for the Plummer sphere. Fig. graphs the 
bound-stars mean radius (short-dash), cluster mean radius (long-dash) and fraction of stars with E > 0 (solid) for the case 
when the SFE = 50%. Comparison with Fig. ^ shows that now the cluster mean radius continues to expand rapidly, while 
the fraction of stars with positive energy and the i? < 0 stellar mean radius flatten out on the same (short) timescale. The 
relatively smooth distribution of stars with F < 0 in these early stages of evolution is a strong indication that the bound stars 
have already settled (mixed) in the inner region (Fig. ^[b,c]). Furthermore, the lack of sharp features in the E > 0 quadrant 
indicates that the motion of the stars involves more orbit-crossing (or, mixing) than for the Plummer model. In other words, 
the assumption that the fraction of bound stars can be identified from (0 independently of the dynamics is now invalidated, 
and we must seek a remedy. 


4.3 An improved scheme 

Our hint comes from Fig. ^ Since the dynamics of the bound stars settles quickly compared with the bulk, we deduce 
that these stars virialise on a very short time scale. We idealise the situation and make this timescale infinitesimal, as for 
the gas-evacuation timescale of the problem. Then a dimensional argument relates the mean square velocity <v^> to the 
self-gravitating mass M and radius R, 


2 \jj lv± 

<V >OC —— OC 0 = 0 * -I- ■ (11) 

H 

Since the escaping stars are removed quickly, and the bound stars evolve on their own timescale, as hinted from Fig. ^ we 
may obtain a better estimate of the bound fraction if we take account of the self-gravitating energy of the bound stars alone 
when defining the velocity DF. Thus let 


^ 4>i,{E < 0 ) 

and truncate the DF at the escape velocity Ue obtained from 

<vl>= 2 X {(t>^,{E < 0) -I- ((igas) . (12) 

The transformation is inexact but useful, since it meets one’s expectation that escaping stars in a Michie-King model will come 
from the outer reach of the cluster (cf. Fig. and hence the transformation (^) corresponds to the limit of a self-gravitating 
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bound core de-coupled from a loose envelope. We may deduce this also from the correlation in energy seen on Fig. @(b,c), 
where stars with initially high-binding energy remained the most bound at later times. 

We repeat our iterative procedure with but with the bounds of integration defined from so that 


1 - Ae = 1 - 


exp [ — 


-erf 4 
2 Vy2 


^ exp [ —- 


1 -3 
- V 
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erf(^ 


- \ 2 


3' 




3/2 


^ 2— :r{(l>-k{E < 0) -f 0gas) ^ 


(13) 


i) = \l < 0) 


The fraction of escapers 1 — Ae is now a stronger function of the computed bound stellar potential 4i*{E < 0). 

The bound fraction computed from (0) is displayed as the solid line on Fig. ^(a) . We find now excellent agreement with 
the nnmerical N-body results. The better fit means we may concentrate on Eq. ( |l3| ) and explore its range of applicability. 
Our intuition in writing down (|^) came from considering the central core as self-gravitating. To see how well the scheme 
performs for less centrally concentrated distributions, we performed a comparison between the analytic scheme and N-body 
calculations of an Michie-King T/cr^ = 3 model. The results are shown on Fig. 1(b)- We have displayed the results obtained 
from (solid line) as well as the old scheme ( 0 ) (dashed line). Once more the N-body data falls much closer to the solid 
line, though the fit is much poorer. The more extended core of the Michie-King T/cr^ = 3 means the region of uniform-density 
stretches closer to the edge of the system: the situation is similar to that of harmonic motion, in which case stars which visit 
the edge escape more easily (Paper I). 

Table ^ lists values of the critical SFE when cluster dissolution would occur, for a sequence of six Michie-King models 
along with their characteristic radii. The parameter T/cr^ ranges from 1 to 12. The bound fraction Ae = M^/M^ is displayed 
on Fig. |^(c), where the free parameter assumes values of 3, 6, 9 and 10. All curves drop to zero bound fraction between 

an SFE e ~ 0.523 for = 3, to « 0.401 for = 9. Models with T/cr^ > 9 reach a higher central density, however the 

mass within the core region is a smaller fraction of the whole; furthermore, the half-mass radius of these models increases with 
T/cr^ so that more stars come close to the system edge: as a result T/cr^ > 9 would dissolve increasingly more easily (Table 


The trend of SEE when dissolution occurs in relation to the core-halo structure of Michie-King models hints that the 
characteristics of the core determine cluster survival. To discover how the mass profile may affect the outcome, we extend our 
iterative scheme to distributions with no harmonic cores. 


5 HERNQUIST PROFILE 

Hernquist (1990) introduced a well-known centrally peaked p oc r~^ density profile. The integrated mass is finite but fills 
the entire space (see the Appendix for details). Experiments with N-body calculations show that a Maxwellian velocity 
distribution, 

f{v)dv^cxv^exp{—v^/2a^)dv, (14) 

provides a self-consistent velocity field everywhere (Hernquist 1993; Boily et al. 2001). Inserting © in yields 


1-Ae(T*) = 1- 


'/T7 exp(-T*) - erf(V2T*) 


(15) 


\/T exp(—T) — erf(\/2T) 

where the dimensionless potential T(r) = <^/cr^(r) and erf(x) is the error function. The relation is once more a function 
of the local stellar potential, T* = . A run of bound fraction Ae decreases with decreasing T*; Ae(T* = 0) = 0 if we 

hold T constant in ©)■ 

To implement the scheme, we truncated the mass profile at a radius r = 200 re, when the spherical volume encloses some 
« 99.5% of the total system mass. Errors in binding energy proved completely negligible for our purpose. The integral ( [l5| ) 
was evaluated as for the King models. Convergence was once more defined as variations in <5(1 — Ae) smaller than one part in 
a million for successive iterations. 

The result of iterating on (^^ is shown as the solid curve on Fig. ^ The curve follows a pattern similar to the solution 
for Plummer or King models but breaks at a lower SFE. Eor example, the fraction of bound stars is still > 50% for an SFE = 
40%. We find that total cluster disruption would take place in this case at a critical SFE = 35.3%, well below the prediction 
(^) and results for Plummer models. 

The velocity dispersion of an Hernquist model drops at the centre of the system as seen when taking the limit a{x —> 0) =0 
in (p^. We would therefore anticipate a smaller fraction of stars escaping from the central region where p oc r~^ compared with 


the enter envelope where p oc r . On Pig. ra(a), we graph the density profiles for Hernquist models and three representative 
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values of the SFE. The equilibrium profile (SFE = 100%) is shown as the solid curve on the figure. The bottom panels graph 
the fraction of bound stars Ae = pt/P* a function of radius. Clearly a greater fraction of stars escapes at large radii than 
at small radii, similarly to Michie-King models (cf. Fig. The fraction of bound stars at the centre Ae(r —> 0) = 0.983 « 1 
for an SFE as low as 37%, while \e{r —> oo) « 0.40 so the ratio > 2 between the two limits. Both Hernquist and Michie-King 
models contrast with the Ae = constant Plummer solution. 


5.1 Other, similar profiles: isochrone and Jaffe models 

Contrary to Michie-King and Plummer models, the Hernquist profile is cuspy in the central region. Is the lower SFE for 
dissolution obtained a result of the central cusp? In order to understand the relative importance of the central cusp, we 
explored the properties of mass profiles that are identical to the Hernquist model at large distances, but differ in the central 
region. We have investigated the properties of two models, the isochrone and Jaffe models, which, like the Hernquist profile, 
have a continuous density profile (no truncation) but finite total mass. Details of each model are given in the Appendix. We 
note here that the density profile of the isochrone equilibrium remains flat in the centre, while the Jaffe model is more steep 
than the Hernsquist model there, pj(r ^ 0) oc r~^. All three show a power-law density p oc at large distances and we 
chose the scales of length so that each model has the same total mass and density at infinity. 

Results are shown for the isochrone (dotted) and the Jaffe (dashed) profiles together with the Hernquist model on Fig. ^ 
It is striking that the SFE when dissolution occurs hardly differs between these three cases, despite different central density 
profiles. We compute for the Jaffe model an SFE = 36.9% and for the isochrone an SFE = 36.4% for total dissolution, which 
compare well with the Hernquist case (35.3%). About 25% of the total mass falls within the central region (inside one length 
scale) in all three cases. This suggests that details of the central mass profile matter little when the mass fraction within 
one core length is comparable in each case. The fact that the DF from which each model derives is a different function of 
energy (see Appendix) also points to a weak dependence of the critical SFE on the shape of the velocity distribution function. 
This point is reinforced if we recall that the Maxwellian distribution (|^) used is not the one derived from Hernquist’s (1990) 
solution. 


5.2 Comparison with Michie-King & Plummer models 

Our experience with Hernquist-type profiles suggests that we should concentrate on averaged values and not dwell on the 
detailed profile of the velocity field or mass. The virial theorem is an obvious reference quantity. Since self-consistent models 
fulfill the scalar virial theorem individually, we would compute for each an absolute ratio of total kinetic to gravitational 
energy = 1/2. Because this holds for equilibrium systems as a whole, this ratio will also be obtained from mean values of 
kinetic and gravitational energy. One way to distinguish between models might be instead to compute a ‘local’ virial ratio, 
taking the average of this ratio over the entire system. We write 


2 \ -1 poo 2 

W\)-mL 


(16) 


where the one-dimensional velocity dispersion a{r) is obtained from integrating the Jeans equation (cf. Eq. ^ below). We 
may compute ( p^ for each one of the models in the previous section and compare with the dissolution SFE. 

The results obtained from ( p^ are listed in Table ^ The quantity {a'^/\(f)\) maps the dissolution SFE for each of the 
Hernquist, isochrone and Jaffe models. For instance, we computed the lowest value {cT/\(l>\) = 0.1805 for the Hernquist model, 
which also allows the lowest SFE before dissolution. The coincidence between the two trends gives us confidence when inter¬ 
preting ((T^/|(()|) as a predictive tool for the critical dissolution SFE. The same approach applied to Michie-King models also 
finds a proportionality between /\(f>\) and the SFE for dissolution (see Table ^). 


The dimensionless ratio ri/(r) provides a diagnostic for the survival of bound clusters in relation to the SFE (Table 

This ratio applied to the family of Michie-King models proved equally reliable as ((t^/|(/|) computed from ( 0 ) to identify 
which model is more robust against mass loss (Table If we extended this conclusion to the three continuous cases of 
Section 5 (Hernquist-type models), we would deduce, wrongly, that the Jaffe model is the most robust of those three : it is in 

fact the least robust (albeit marginally so), while ri/{r) = 0.169 is easily the smallest value for these models (Table ^ rggg 

2 

is the radius enclosing 99.9% of the mass). Hence the mass profile alone, and by implication the potential, does not provide a 
criterion for survival applicable to a wide variety of cases. The average ratio (|l6|) appears a better diagnostic, at least as far 
as we are concerned with models similar in character to one another. 


The case of the Plummer sphere cautions against applying ( 0 ) to arbitrary systems. In this case (ct^/|<(>|) = 1/6 is lower 
than that computed for the isochrone, Hernquist or Jaffe models, yet its dissolution SFE is larger than that of any of these 
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three models. It is worth pointing out that at large SFE, the Plummer sphere actually is more robust than any of the other 
models considered (see Fig. ^[c,d]). The swift transition to low M*/M* bound fraction and eventual dissolution obtained for 
the Plummer may instead indicate peculiarities in the dynamics of expansion for this model: we noted that the central volume 
of Michie-King models re-collapses on itself on a short timescale, while the Plummer solution does not (Figs. and 


6 TIDAL AND COLLISIONAL EFFECTS 

General comments are possible based on data listed in Tables ^ to These list the factors by which concentric shells expanded 
in the N-body experiments carried out. Clearly the equilibrium clusters that ensued suffered much expansion and mass loss. 
The final equilibrium mass is found from integrating 


M* = / F{r,v) d^r d^w oc (/)r*(uj)^^^ 

Jv^, 

where (/) is averaged over the entire phase-space volume V* and 


2 

oc 


GM^ 

r* 


since the cluster achieved a new virial equilibrium by hypothesis. Isolating for (/) we find 


(/) (X . 

From e.g. Table ^ we find for an SFE = 50% a bound stellar mass of ~ 31% of the total initial mass, while the radius now 
expanded on the mean by a factor ~ 5, relatively large compared with, e.g., Scorpius data, which suggest an expansion factor 
of a few between clusters of different ages (Brown 2001). Inserting this in the relation above, we find a net decrease in the mean 
phase-space density by a factor \/3/5®^^ ~ 6“^. The rough derivation shows that although a star cluster may survive gas loss, 
its binding energy will be much lowered and therefore the cluster will be subject to significant heating from an external tidal 
potential. None of the DF’s considered here account for this heating. This is surely at work in the case of open clusters in the 
Galaxy. Gare must be taken when making predictions for the short-term evolution of clusters filling their Roche tidal radius 
however (Terlevich 1987). Baumgardt (2001) finds the orbits of E « 0 tidally unbound stars persisting for several cluster 
dynamical times, owing to the shared galactic orbital motion by the star and the cluster centre of mass. Thus unbound stars 
seen in projection against the sky may still be counted as cluster members and contribute to the potential energy, while the 
energy criteria we have used assumes instant removal of all stars with E > 0. This would offset to a certain degree the net 
effect of tides due to the large expansion factors seen in the numerical computations. Overall, tidal fields will strip clusters of 
more stars than calculated from (^ , however the quantities may only be worked out once cluster orbital parmeters are specified. 


The dynamics of Michie-King models following gas expulsion proceeds on a short timescale in their centre. We found 
better agreement with numerical experiments when treating the stars with E < 0 initially as a self-gravitating group, so re¬ 
setting the cut-off velocity of the DF in the problem posed. By implications, the central region undergoes much shell crossing 
and hence collisional effects there may yet be of importance in this problem. Kroupa, Aarseth & Hurley (2001) evolved an 
embedded Plummer model with a high-precision direct-summation A-body code with delayed, but then near-to instantaneous, 
gas-removal. They find a fraction ~ 30% of stars remain bound despite a low SFE of « 33%(e = 0.33). The results are similar 
to the analytical results for the Hernquist model (Eig. solid curve), while we would have expected complete dissolution 
for strictly collision-less evolution for this SFE. The delay in removing the gas allows the central part to increase its binding 
energy by two-body relaxation. The central part therefore increases its binding energy at the expense of the outer envelope. 
The effect reported there, combined with the result of the present work, help understand the conditions required to form 
bound objects when the SFE ~ 20% or lower (Lada 1999). 


7 GENERAL CONCLUSIONS 

We carried out numerical N-body calculations to verify that the iterative scheme introduced in Paper I does not grossly under¬ 
estimate the effects of dynamical evolution: We found excellent agreement between numerical calculations and analysis for the 
case of Plummer spheres. Rapid virialisation of the central region of Michie-King models required a more careful selection of the 
bound stars from the velocity DF for agreement (cf. Eq. ^). This confirms our conclusion of Paper I, namely that the initial 
DF of embedded clusters defines the stellar population of bound, gas-free clusters. Some highlights from the present study are: 
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1) The mass-weighted mean ratio ((t^/|(^|) equation ( p^ provides a diagnostic to predict critical SFE for total disruption. 
Thus systems with a lower value computed from will dissolve at a lower SFE. We have found this to match the analytic 
results for all DF considered with the exception of the Plummer DF (Table ^). 

2) When a self-gravitating system evolves such that individual star-star encounters are ignored, any member of the 
Michie-King family of star cluster profile will disrupt if less than 40% of the gas is turned into stars initially. 

3) A peaked Hernquist-type prohle with decaying central velocity dispersion provides the most robust configuration 
against gas loss amongst the DF we considered. We found similar values of critical SFE for Hernquist, Jaffe and isochrone 
DF, which all are predicted to dissolve when the SFE is less than « 35%. 

The results obtained here gives a clue to understanding the demographics of massive stars in cluster cores, where the 
more massive stars are found preferentially in denser cores (Testi et al. 1997, 1999): indeed only the more cencentrated clusters 
would survive rapid gas expulsion driven by high-mass stars, while others would disperse. The presence of a population of 
expanding stars around such dense cores would be a strong indication that the clusters underwent such a phase of gas removal. 
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APPENDIX: DETAILS OF THE DF’S 

In this appendix we give details of the numerical integrations for the DF’s of the three models with mass densities p{r) 
satisfying p(r ^ oo) oc discussed in Section 
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Hernquist DF 

The density p and potential of this model vary radially according to 


p{r or x) = 


M 


M 1 


1 




( 17 ) 


2-k r (r + Tc)® 2nr^ x {x + 1)® 2nGr'^ x{x + 1)^ ’ 
with Tc a free length fixing the point of the power-law turnover, and x = rjrc- 

The velocity field may be recovered in two ways. A useful approximation based on the first moment of Jeans equation 
(cf. BT-l-87, §4.2; Hernquist 1993) is to evaluate the one-dimensional isotropic velocity dispersion from 


1 7°° 

^(r) = y p{x)V(l>{x) dx 


which may be computed from (|17D; we find 


o-^{x) = 4>{x) a; (1 -I- x)‘'‘ < In 


1 -I- a; 


1/4 


1/3 


{l + xY (l-|-a;)3 


1/2 1 1 
(1 -I- xY 1 + X j ’ 


(18) 


(19) 


which constrains f{v) locally through a{r or a:); we thus have the freedom to choose any profile satisfying a(r). A Maxwellian 
velocity DF is found to give stable equilibria in N-body realisations of the model when the DF is cut off at the escape velocity 
Ve, and a re-normalised to (Hernquist 1993) 


cr^(r) = <T^(r)- 


^erf(ne/t/2cr) - -^^exp{-v'i/2a^) 


|-erf(ue/\/2cr)- exp{—Ve/2a^) -— exp{—v^/2a) 


2^f2o 


6\/2o 


such that (pY) — 3(t^ as for a Maxwellian profile. 

The velocity field may be recovered directly from the model DF. The DF which gives rise to (j^) is obtained directly 
from an Abel transformation (Hernquist 1990) 


8^7r3 


M 1 I 3 sin ^ q + q\/^- <1^(1 “ 2g^)(8g'‘ - - 3) 


2x5/2 

(1-g ) 


where the definitions 


2 GM 
Vg = - ; q = 

Tc 


Exc 

'Wid 


= U>- 


1/2 


; l>(r) = ||<E>(r) 


apply. The bound stellar fraction Ae is then evaluated from 


Ae = 


pu 

/ F[E]t 
Jo 


1x^2 l-'*/(l>-l) 
where we set u = \/ 2e$ < \/ 2$. The indefinite integral then reads 


p U 

/ F[E] 
Jo 


V dv = 


sm-\^l^^^) + (^-2 + 2^-u^^ y 6<l'‘ + ll$^u^ (2-b u^) - l>® (6-bl7u^)-b 

+ {l(i-l) (2-2* + . 


lY (—3 -b -b 2 u"^) — 2 $ (n^ -b 6 n"* -b u®) 

The solution Xe{r) is obtained by repeated iterations, upon a re-evaluation of the potential. The upper bound of integration 
is then u = xj 23>*(r), i.e. root of the dimensionless escape velocity derived from the stellar potential. 


Jaffe DF 

The model discussed by Jaffe (1983) derives from the potential-density pair 
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pj(r) = -- • $j(r) = In 

47rr? r'^{r + rj)'^ ’ Vj 


r + r 


In the above expression rj is a free scale length. Near the centre we find pj{r ^ 0) oc r ^ is diverging like a power of r. At 
large distances 

hm pj(r) = -— — 

r—*oo 47r 

which is identical to the density of an Hernquist model in the same limit if we make rj = 2rc. Both systems then have the 
same total integrated mass, M. The velocity dispersion o-(r) is known explicitly everywhere. Let x = rjrj. Then 




aj{x) = -{1 -4a:- 18a:^ - 12a;®} - 12a:®(l + a:®)^<E>(a;) 


constant. 


Thus the velocity dispersion levels off near the centre. The DF may be expressed explicitly for this model (Jaffe 1983; BT+87) 

M 


Fj{E) = 


X {f_(v^) - V2F.{^) - V2F+{^) + F+{y/2q)] 


7r3 

with q defined as before but with the substitution —> Xj. The dimensionless integral F±{u) is a modified Dawson integral 


^^( m )=f dy. 

Jo 

The substitution Fj{F) above in the integral (^) yields no simple expression and we must proceed numerically. The integrals 
F± however are evaluated efficiently with numerical integrals (cf. Press et al. 1992, §6.10). 


Isochrone DF 

The potential-density pair of the isochrone model are (BT-l-87) 


with the definition a{r) 


^ 3a®(r)(b-ha(r))-r®(b-ha(r)) . * . x GM 

^ 47ra3(r)(6-f a(r))3 ’ ^ 6-|-a(r) 


62 -h 


r2. The asymptotic limit pi{r —> 0) oc constant, whereas at large distances 


lim pi{r) = 


We therefore recover the same density at large distances for all models if their total integrated mass M is equal and the 
lengths b = Xc = d /2. In terms of the parameter q the DF is known explicitly from (Henon 1960) 


Unfortunately we found no algebraic expression upon inserting Fi in (|^ and integrating. Numerical integration presented no 
difficulty. 
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Table 1. Computed quantities for Plummer models together with results of N-body calculations. 


SFE Q Bound fraction Ri,/Ro Ratio of Lagrangian radii 


100 X e 


Eq# 

Numerical 

Eq (1) 

Numerical 

10% 

20% 

30% 

50% 

70% 




±0.05 


(mean) 






100% 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

90% 

1.11 

0.99 

1.00 

1.13 







80% 

1.25 

0.99 

0.98 

1.33 

1.37 

1.47 

1.38 

1.33 

1.30 

1.35 

70% 

1.43 

0.99 

0.95 

1.75 

1.68 

1.74 

1.64 

1.60 

1.62 

1.80 

60% 

1.67 

0.95 

0.86 

3.00 

2.40 

2.21 

2.11 

2.10 

2.44 

3.12 

55% 

1.82 

0.92 

0.78 

5.50 

3.39 

2.69 

2.61 

2.71 

3.68 

5.26 

50% 

2.00 

0.85 

0.66 

OO 

4.35 

3.29 

3.38 

4.13 

6.63 

n.a. 

45% 

2.22 

0.69 

0.44 


11.5 

5.50 

10.6 

18.4 

n.a. 

n.a. 

40% 

2.50 

0.00 

0.00 


OO 







Table 2. Computed quantities for a King = 9 model together with results of N-body calculations, n.a. = non-available. 


SFE Q Bound fraction R*/i?o Ratio of Lagrangian radii 


100 X e 


Eq® 

Numerical 

Eq (|) 

Numerical 

10% 

20% 

30% 

50% 

70% 



±0.05 

(mean) 






100% 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

90% 

1.11 

0.99 

1.00 

1.13 







80% 

1.25 

0.97 

0.98 

1.33 

1.35 

1.45 

1.36 

1.31 

1.29 

1.33 

70% 

1.43 

0.91 

0.92 

1.75 

1.65 

1.66 

1.64 

1.58 

1.63 

1.75 

60% 

1.67 

0.79 

0.80 

3.00 

2.31 

1.96 

2.10 

2.15 

2.34 

3.00 

50% 

2.00 

0.63 

0.64 

OO 

5.56 

3.12 

3.74 

4.89 

10.5 

n.a. 

40% 

2.50 

0.27 

0.39 


7.18 

4.65 

6.79 

10.1 

n.a. 

n.a. 

30% 

3.33 

0.00 

(0.10) 


OO 






20% 

5.00 

0.00 

0.00 









Table 3. Computed quantities for a King = 3 model together with results of N-body calculations, n.a. = non-available. 


SFE Q Bound fraction Ri,/Ro Ratio of Lagrangian radii 


100 X e 


Eq (y) 

Numerical 

Eq (1) 

Numerical 

10% 

20% 

30% 

50% 

70% 




±0.05 


(mean) 






100% 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 

90% 

1.11 

0.99 

1.00 

1.13 







80% 

1.25 

0.94 

0.99 

1.33 

1.25 

1.26 

1.26 

1.25 

1.25 

1.25 

70% 

1.43 

0.86 

0.97 

1.75 

1.55 

1.53 

1.52 

1.52 

1.54 

1.63 

60% 

1.67 

0.68 

0.84 

3.00 

2.18 

1.93 

1.95 

1.98 

2.13 

2.93 

50% 

2.00 

0.00 

0.53 

OO 

3.09 

2.92 

3.05 

3.30 

n.a. 

n.a. 

40% 

2.50 

0.00 

(0.11) 


> 19.00 

> 19.00 

n.a. 

n.a. 

n.a. 

n.a. 

30% 

3.33 

0.00 

0.00 


OO 
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Table 4. Comparison of model characteristics. The SFE at which dissolution occurs is given for each model along with mean and 
half-mass lengths ri for each case. The radius rggg within which the mean radius was evaluated defines the sphere enclosing 99.9% of 

the mass. The characteristic lengths I are given for each case. For Michie-King models the core length Vq defined by (j^ is set such that 
the cluster radius rt =5/4. 


King 

length 1 

<r> 

n 

r 1 

r 1 / <r> 


SFE 

= 

To = 



2 

2 

Eq. 

e 

1 

0.58 

0.463 

1.250 

0.440 

0.950 

0.808 

0.529 

3 

0.20 

0.377 

1.250 

0.343 

0.911 

0.786 

0.523 

6 

0.10 

0.261 

1.250 

0.203 

0.777 

0.709 

0.487 

9 

2.50 X 10-3 

0.316 

1.250 

0.229 

0.725 

0.685 

0.401 

10 

1.25 X 10-3 

0.379 

1.250 

0.303 

0.801 

0.722 

0.449 

12 

4.98 X lO-"' 

0.479 

1.250 

0.426 

0.893 

0.764 

0.496 

Model 

length 1 

<r> 

^999 

r 1 

r 1 / <r> 


SFE 





2 

2 

Eq. 

e 

Isochrone 

6 = 1 

12.83 

1999.3 

3.060 

0.238 

0.245 

0.364 

Hernquist 

rc = l 

12.22 

1998.5 

2.414 

0.198 

0.181 

0.353 

Jaffe 

rj = 2 

11.83 

1998.0 

2.000 

0.169 

0.252 

0.369 

Plummer 

Rp = 1 

1.925 

38.714 

1.305 

0.678 

0.1667 

0.442 



time [m.u.] 



E[t = 0] 



E[t = 0] 


Figure 1. Time-evolution of an N-body calculation. The model Plummer cluster had parameters Q = 2 and a truncation radius 
R = 22 Rp = 1.25 m.u.; N = 50, 000 particles were used in the calculation, (a) This graphs the mean radius of all stars (long-dash) and 
stars with E < 0 (short dash) as function of time. (One m.u. of time ^ 1.3 £cr defined in [Ml.) The solid line shows the fraction of stars 
N{E > 0)/^” with positive energy, (b) and (c): particle energy versus their initial energy (after gas expulsion) at times £ = 1/2 and 2, 
respectively. The triangle indicates stars that had negative energy initially which are now unbound. 
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r / r 
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Figure 2. Initial and final (after iterations) density profiles for (a) Plummer and (b) King ^ jcP' = 9 models for three values of the SFE, 
e. The run of Ae, i.e. the ratio of density to the initial density, is also shown as a function of radius (bottom panels). Note how the King 
profile becomes steeper for small e, while a Plummer model loses mass equally at all radii. 
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s.f.e. s.f.e. 




s.f.e. s.f.e. 


Figure 3. Ratio Ae = M^jM* of bound to initial stars as function of star formation efficiency e. (a) The results for King ^/cj^ = 9 
models derived from (^) are shown as the dashed line; the improved solution ( [l^ ) is the solid line. The results of numerical N-body 
computations are shown as open circles and listed in Tables ^ and ^ (b) As (a), but for =3 King models, (c) Solutions for four 

different King models with W = ^ = 3, 6, 9 and 10. (d) Comparison between two King models and a Plummer model shown as solid. 

The Plummer model dissolves for an SFE below ^ 0.44. 
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time [m.u.] E [ t = 0 ] 

Figure 4. Time-evolution of an N-body calculation. As for Fig. 1 but for a Michie-King model with parameters ^ = 9 and Q 

(One model unit of time 1.3 ter defined in 0.) 
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s.f.e. 


Figure 5. Bound fraction M*/M* versus SFE for three cuspy models (no harmonic cores) identified in the key. All three curves break 
sharply below SFE = 0.40. 
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r / r 


c 


Xe 1 
0.8 
0.6 
0.4 
0.2 


0.01 0.10 1.00 10.00 




0.01 


0.10 


1.00 


10.00 


Figure 6. Density profiles (top) and fraction of bound stars (bottom) for two models, (a) The Hernquist profile peaks at the centre, 
p (X whereas on (b) the isochrone model has p(r ^ 0) = constant. Both have p oc at large distances. When comparing the two 
models, we find the bound fraction Ae assumes similar values at large distance for the same SFE; the bound fraction remains ^ 1 at the 
centre even for low SFE in the case of the Hernquist model, but drops significantly for the isochrone solution. 



















